############################################
###REPLICATION CODE FOR ASH (2023)##########
###Fallow agriculture and climatic stress###
###independently predict migration during###
###Syria's 2006-10 Drought##################
###forthcoming in Regional Environmental ###
###Change###################################

library(foreign)
library(readstata13)
library(list)
library(tidyverse)

library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #still crashes
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(kfigr) #figure referencing for markdown
library(knitr) #knitting
library(pander) #pandering
library(geosphere) #distance tools
library(Rcpp) #C++ tools
library(RcppArmadillo) #C++ tools
library(viridis) # color map
library(multiwayvcov) # multiway clustered standard errors for lm
library(gtable) # for switching facet labels in ggplot
library(matrixStats) # for rowMins function
library(readstata13) # read in newer stata files
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(modelr)

setwd("C:/Users/kash8/OneDrive/Documents/GIS/Eklund et al 2022 Extension")
setwd("E:/OneDrive/Documents/GIS/Eklund et al 2022 Extension")

###setwd("~Ash 2023 REC Replication")

master_fallow<-read.dta("master_fallow.dta")
master_fallow$name_en<-as.factor(master_fallow$name_en)
master_fallow$admin1<-as.factor(master_fallow$admin1)
master_fallow$admin2<-as.factor(master_fallow$admin2)
master_fallow <- master_fallow[!grepl("Arbin|Damascus|Hajar Aswad|Jaramana", master_fallow$name_en), ]

###Figures 1a-c: FE OLS of Year-of, 2008 and 2000 Fallowness on Year-of nighttime lights####

###Annual####
drought.2000 <- felm(lights ~ fallow | 0 | 0 | admin1, data=master_fallow[master_fallow$year==2000,])
summary(drought.2000)
drought.2001 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2001,])
summary(drought.2001)
drought.2002 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2002,])
summary(drought.2002)
drought.2003 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2003,])
summary(drought.2003)
drought.2004 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2004,])
summary(drought.2004)
drought.2005 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2005,])
summary(drought.2005)
drought.2006 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2006,])
summary(drought.2006)
drought.2007 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2007,])
summary(drought.2007)
drought.2008 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2008,])
summary(drought.2008)
drought.2009 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2009,])
summary(drought.2009)
drought.2010 <- felm(lights ~ fallow| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2010,])
summary(drought.2010)

###2008 as variable
drought.2009.08 <- felm(lights ~ fallow_2008| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2009,])
summary(drought.2009.08)
drought.2010.08 <- felm(lights ~ fallow_2008| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2010,])
summary(drought.2010.08)
###2000 as variable
drought.2001.00 <- felm(lights ~ fallow_2000| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2001,])
summary(drought.2001.00)
drought.2002.00 <- felm(lights ~ fallow_2000| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2002,])
summary(drought.2002.00)
drought.2003.00 <- felm(lights ~ fallow_2000| 0 | 0 | admin1, data=master_fallow[master_fallow$year==2003,])
summary(drought.2003.00)
drought.2004.00 <- felm(lights ~ fallow_2000| 0 | 0 | admin2, data=master_fallow[master_fallow$year==2004,])
summary(drought.2004.00)
drought.2005.00 <- felm(lights ~ fallow_2000| 0 | 0 | admin2, data=master_fallow[master_fallow$year==2005,])
summary(drought.2005.00)


fallow_year<- read.csv("fallow_year.csv",sep = "\t")


FELM <- data.frame(Treatment = fallow_year$year,
                         Difference = fallow_year$coef,
                         SE = fallow_year$se,
                         modelName = "Coefficient") 

FullModelFrame <- data.frame(rbind(FELM))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(y = Treatment, xmin = Difference - SE*interval1,
                                xmax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(y = Treatment, x = Difference, xmin = Difference - SE*interval2,
                                 xmax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("\u0394 in Nahiya Nighttime Lights") + ylab("")+ xlim(-20,20)+ theme(legend.position="none") #+ labs(color='Model')
#scale_y_discrete(labels = function(y) str_wrap(y, width = 10))+
#  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
#       panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 + scale_color_manual(values=c("#000000"))+scale_y_reverse(breaks=c(2010,2009,2008,2007,2006,2005,2004,2003,2002,2001,2000))
print(zp1) # The trick to these is position_dodge().

fallow_2008<- read.csv("fallow_2008.csv",sep = "\t")


FELM2008 <- data.frame(Treatment = fallow_2008$year,
                   Difference = fallow_2008$coef,
                   SE = fallow_2008$se,
                   modelName = "Coefficient") 

FullModelFrame <- data.frame(rbind(FELM2008))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(y = Treatment, xmin = Difference - SE*interval1,
                                xmax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(y = Treatment, x = Difference, xmin = Difference - SE*interval2,
                                 xmax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("\u0394 in Nahiya Nighttime Lights") + ylab("")+ xlim(-15,20)+theme(legend.position="none",axis.text.y=element_blank()) #+ labs(color='Model')
#scale_y_discrete(labels = function(y) str_wrap(y, width = 10))+
#  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
#       panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 + scale_color_manual(values=c("#000000"))+scale_y_reverse(breaks=c(2010,2009,2008,2007,2006,2005,2004,2003,2002,2001,2000))
print(zp1) # The trick to these is position_dodge().

fallow_2000<- read.csv("fallow_2000.csv",sep = "\t")


FELM2000 <- data.frame(Treatment = fallow_2000$year,
                   Difference = fallow_2000$coef,
                   SE = fallow_2000$se,
                   modelName = "Coefficient") 

FullModelFrame <- data.frame(rbind(FELM2000))
#model1Frame$Treatment <- factor(model1Frame$Difference, as.character(model1Frame$Difference))
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(FullModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) 
zp1 <- zp1 + geom_linerange(aes(y = Treatment, xmin = Difference - SE*interval1,
                                xmax = Difference + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(y = Treatment, x = Difference, xmin = Difference - SE*interval2,
                                 xmax = Difference + SE*interval2),
                             lwd =1, position = position_dodge(width = 1/2))
zp1 <- zp1  + theme_bw(base_size = 14) + xlab("\u0394 in Nahiya Nighttime Lights") + ylab("")+ xlim(-15,20)+theme(legend.position="none",axis.text.y=element_blank()) #+ labs(color='Model')
#scale_y_discrete(labels = function(y) str_wrap(y, width = 10))+
#  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
#       panel.background = element_blank(), axis.line = element_line(colour = "black"))
zp1 <- zp1 + scale_color_manual(values=c("#000000"))+scale_y_reverse(breaks=c(2010,2009,2008,2007,2006,2005,2004,2003,2002,2001,2000))
print(zp1) # The trick to these is position_dodge().


##################################################################################
####Figures 1d and 1e: Maps of fallowness and 2007-10 Lights Change###############
##################################################################################

library(gstat)
library(tmap)


fallow2008<-readOGR(dsn=".",layer="SYR_ADM3_agriculture")
lights<-readOGR(dsn=".",layer="SYR_ADM3_lights")

syriabase<-readOGR(".","SYR_adm0")
plot(syriabase,col="lightgray")

fallow2008@data[["MEAN"]]
fallow2008@data$MEAN[fallow2008@data$MEAN ==0] <-NA
tm_shape(syriabase)+tm_fill()+tm_shape(fallow2008)+ tm_borders()+tm_fill(col = 'MEAN',title = "Fallowness",n=10, palette="Oranges")+tm_layout(legend.outside=TRUE,legend.outside.size = .15,frame = FALSE)
tm_shape(syriabase)+tm_fill()+tm_shape(lights)+ tm_borders()+tm_fill(col = 'MEAN_12',title = "Nighttime Lights",n=10, palette=c("black","white"))+tm_layout(legend.outside=TRUE,legend.outside.size = .15,frame = FALSE)
lights@data$lights0710<-lights@data$MEAN_12_13-lights@data$MEAN_1
tm_shape(syriabase)+tm_fill()+tm_shape(lights)+ tm_borders()+tm_fill(col = 'lights0710',title = "Light Change",n=10, palette=c("black","white"))+tm_layout(legend.outside=TRUE,legend.outside.size = .15,frame = FALSE)


###Calculating Nighttime lights Saturation with 2010 DMSP-OLS Data#######

nightlights2010 <- raster('F182010.v4d_web.stable_lights.avg_vis.tif')
syr_adm3_noflares2<-readOGR(dsn=".", layer="SYR_adm3_noflares")
plot(syr_adm3_noflares2)
subset_fallow<-master[ which(master$year==2010),]
subset_fallow<-rename(subset_fallow, name = NAME_EN)
subset_fallow<-rename(subset_fallow, admin1 = ADMIN1_EN)
subset_fallow<-rename(subset_fallow, admin2 = ADMIN2_EN)

subset_fallow<-within(subset_fallow, rm(year,X,temp,prcp,mean.temp,mean.prcp,lights,fallow_2000,fallow_2008))
syr_adm3_noflares2@data<-merge(syr_adm3_noflares2@data,subset_fallow,all.x=TRUE,allow.cartesian=TRUE)
syr_adm3_noflares3<-syr_adm3_noflares2[!is.na(syr_adm3_noflares2@data$fallow) ,]
plot(syr_adm3_noflares3)

fallow_lights2010 <- crop(nightlights2010, extent(syr_adm3_noflares3))
fallow_lights2010 <- mask(fallow_lights2010, syr_adm3_noflares3)

freq(fallow_lights2010,useNA='no')

#1,377 pixels are 60 or over from a total of 163,630 pixels --  



#####Table 1: Temperature, Precipitation and Fallowness in Time-by-Location model on Nighttime Lights#####


master<-read.csv("master_rev2.csv", sep = ';')
master_merge<-read.dta("master_merge.dta")
master<-merge(master,master_merge,allow.cartesian=TRUE)
master_merge$name<-as.factor(master_merge$name)
master_merge <- master_merge[!grepl("Arbin|Damascus|Hajar Aswad|Jaramana", master_merge$name), ]

# year prior deviations----
setDT(master)
setkey(master, name, year)
master[, prcp.dev := shift(prcp - mean.prcp), by=.(name)]
master[, temp.dev := shift(temp - mean.temp), by=.(name)]


drought.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow | name  + as.factor(admin1):as.factor(year) | 0 | admin1, data=master[master$year >= 2006 & master$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow, data=master[master$year >= 2006 & master$year <= 2010,])
predrought.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow | name  + as.factor(admin1):as.factor(year) | 0 | admin1, data=master[master$year < 2006,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow, data=master[master$year < 2006,])
all.felm_full <- felm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow | name + as.factor(admin1):as.factor(year) | 0 | admin1, data=master, exactDOF=TRUE, psdef=FALSE)
all.lm_full <- lm(lights ~ prcp.dev + temp.dev + prcp.dev:temp.dev+fallow, data=master)

summary(drought.felm_full)
summary(drought.lm_full)
summary(predrought.felm_full)
summary(predrought.lm_full)
summary(all.felm_full)
summary(all.lm_full)

###Plot Table 1###

stargazer(all.felm_full, predrought.felm_full, drought.felm_full,
          type='latex',
          title = 'Climatic Stress and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', 'Fallow Agr. Pct.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are multiway clustered on','second-level administrative districts and years.',
                    'Only the 185 nawahi Eklund et. al. 2022', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)




##########################
####APPENDIX##############
##########################

###Models for Effect of Fallowness on Nighttime Lights
drought.fallow.felm <- felm(lights ~ fallow | name_en + as.factor(admin1):as.factor(year) | 0 | admin2, data=master_fallow[master_fallow$year >= 2008  & master_fallow$year <= 2010,], exactDOF=TRUE)
drought.fallow.somefe <- felm(lights ~ fallow | admin1 + year | 0 | admin2, data=master_fallow[master_fallow$year >= 2008 & master_fallow$year <= 2010,])
predrought.fallow.felm <- felm(lights ~ fallow | name_en + as.factor(admin1):as.factor(year) | 0 | admin2, data=master_fallow[master_fallow$year < 2008,], exactDOF=TRUE)
predrought.fallow.somefe <- felm(lights ~ fallow | admin1 + year | 0 | admin2, data=master_fallow[master_fallow$year < 2008,])
all.fallow.felm <- felm(lights ~ fallow | name_en + as.factor(admin1):as.factor(year) | 0 | admin2, data=master_fallow, exactDOF=TRUE)
all.fallow.somefe <- felm(lights ~ fallow | admin1 + year | 0 | admin2, data=master_fallow)

summary(drought.fallow.felm)
summary(drought.fallow.somefe)
summary(predrought.fallow.felm)
summary(predrought.fallow.somefe)
summary(all.fallow.felm)
summary(all.fallow.somefe)



###Table A1: Limited fixed effects for Fallowness on Nighttime Lights ###

stargazer(all.fallow.somefe, predrought.fallow.somefe, drought.fallow.somefe,
          type='latex',
          title = 'Agricultural Fallowness and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c( 'Fallow Agr. Pct.'),
          add.lines = list(c("First-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "No", "No", "No")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2022)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)
####Table A2: Time-by-location Fixed Effects Models for Fallowness on Nighttime Lights####

stargazer(all.fallow.felm, predrought.fallow.felm, drought.fallow.felm,
          type='latex',
          title = 'Agricultural Fallowness and Nighttime Lights',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c( 'Fallow Agr. Pct.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2022)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)


#####Table A3: Lagged Fallowness on Nighttime Lights Time-by-Location####

setDT(master)
setkey(master, name, year)
master[, fallow_lag := shift(fallow), by=.(name)]
master[, fallow_future := shift(fallow, type='lead'), by=.(name)]

drought.fallow.felm.lag<- felm(lights ~ fallow+fallow_lag | name + as.factor(admin1):as.factor(year) | 0 | admin2, data=master[master$year >= 2008  & master$year <= 2010,], exactDOF=TRUE)
predrought.fallow.felm.lag <- felm(lights ~ fallow+fallow_lag | name + as.factor(admin1):as.factor(year) | 0 | admin2, data=master[master$year < 2008,], exactDOF=TRUE)
all.fallow.felm.lag  <- felm(lights~ fallow+fallow_lag | name + as.factor(admin1):as.factor(year) | 0 | admin2, data=master, exactDOF=TRUE)

summary(drought.fallow.felm.lag)
summary(predrought.fallow.felm.lag)
summary(all.fallow.felm.lag)

stargazer(all.fallow.felm.lag, predrought.fallow.felm.lag, drought.fallow.felm.lag,
          type='latex',
          title = 'Agricultural Fallowness and Nighttime Lights with Lag',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Nighttime Lights',
          covariate.labels = c( 'Fallow Agr. Pct.', 'Fallow Agr. Pct. Lagged 1 Year'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2022)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

#####Table A4: Effect of Lagged Nighttime Lights on Fallowness####


drought.fallow.felm.rev <- felm(fallow_future ~ lights | name  + as.factor(admin1):as.factor(year) | 0 | admin2, data=master[master$year >= 2007  & master$year <= 2010,], exactDOF=TRUE)
predrought.fallow.felm.rev <- felm(fallow_future ~ lights | name  + as.factor(admin1):as.factor(year) | 0 | admin2, data=master[master$year < 2007,], exactDOF=TRUE)
all.fallow.felm.rev  <- felm(fallow_future~ lights | name + as.factor(admin1):as.factor(year) | 0 | admin2, data=master, exactDOF=TRUE)

summary(drought.fallow.felm.rev)
summary(predrought.fallow.felm.rev)
summary(all.fallow.felm.rev)

stargazer(all.fallow.felm.rev, predrought.fallow.felm.rev, drought.fallow.felm.rev,
          type='latex',
          title = 'Nighttime Lights and Future Agricultural Fallowness',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Fallowness in Following Year',
          covariate.labels = c( 'Nighttime Lights'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are clustered on','second-level administrative districts.',
                    'Only the 185 nawahi Elklund et. al. (2022)', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="30pt"
)






###Table A5: Time by location fixed effects of temperature and precipitation on fallow agriculture#####

drought.felm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name  + as.factor(admin1):as.factor(year) | 0 | admin1, data=master[master$year >= 2006 & master$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
drought.lm_fallow <- lm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=master[master$year >= 2006 & master$year <= 2010,])
predrought.felm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name + as.factor(admin1):as.factor(year) | 0 | admin1, data=master[master$year < 2006,], exactDOF=TRUE, psdef=FALSE)
predrought.lm_fallow <- lm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=master[master$year < 2006,])
all.felm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name  + as.factor(admin1):as.factor(year) | 0 | admin1, data=master, exactDOF=TRUE, psdef=FALSE)
all.lm_fallow <- lm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev, data=master)

summary(drought.felm_fallow)
summary(drought.lm_fallow)
summary(predrought.felm_fallow)
summary(predrought.lm_fallow)
summary(all.felm_fallow)
summary(all.lm_fallow)

###Display Table A5####

stargazer(all.felm_fallow, predrought.felm_fallow, drought.felm_fallow,
          type='latex',
          title = 'Climatic Stress and Fallow Agriculture',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Agricultural Fallowness',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "Yes", "Yes", "Yes")),
          notes = c('Standard errors are multiway clustered on','second-level administrative districts and years.',
                    'Only the 185 nawahi Eklund et. al. 2022', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

####Table A6: Third Level Admin and Year Fixed Effects of Lagged Drought on Fallowness####



drought.somefelm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name + year  | 0 | admin1, data=master[master$year >= 2006 & master$year <= 2010,], exactDOF=TRUE, psdef=FALSE)
predrought.somefelm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name + year | 0 | admin1, data=master[master$year < 2006,], exactDOF=TRUE, psdef=FALSE)
all.somefelm_fallow <- felm(fallow ~ prcp.dev + temp.dev + prcp.dev:temp.dev | name + year | 0 | admin1, data=master, exactDOF=TRUE, psdef=FALSE)

summary(drought.somefelm_fallow)
summary(predrought.somefelm_fallow)
summary(all.somefelm_fallow)


stargazer(all.somefelm_fallow, predrought.somefelm_fallow, drought.somefelm_fallow,
          type='latex',
          title = 'Climatic Stress and Fallow Agriculture',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Agricultural Fallowness',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes"),
                           c("First-Level Admin:Year FE", "No", "No", "No")),
          notes = c('Standard errors are multiway clustered on','second-level administrative districts and years.',
                    'Only the 185 nawahi Eklund et. al. 2022', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)


###Figure A1:Visualizing Interaction Effect of Lagged Precipitation and Temperature on Fallowness; Limited FEs####

continuousinteraction <- function(felmmodel, lmmodel, rhs1, intxn, minimum="min", maximum="max", miny=-8, maxy=1, incr="default", num_points = 1000, conf=.95, title="Interaction marginal effects plot", xlabel="Value of Intxn", ylabel="Estimated Marginal Effect", maincolor='black', secondcolor='black'){
  
  # Get coefficients
  beta_1 = felmmodel$coefficients[row.names(felmmodel$coefficients)==rhs1]
  
  # Get coefficient if interaction term==100
  beta_3 = felmmodel$coefficients[row.names(felmmodel$coefficients)==paste0(rhs1, ':', intxn)]
  
  # Extract Variance Covariance matrix
  covMat = felmmodel$clustervcv
  
  # Extract the data frame of the model
  mod_frame = model.frame(lmmodel)
  
  # Set range of the interaction variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[intxn]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[intxn]])
  }else{
    max_val = maximum
  }
  
  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum X value greater than maximum value.")
  }
  
  # Determine intervals between values of the x variable
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }
  
  # Create list of intxn values at which marginal effect is evaluated
  intxn_2 <- seq(from=min_val, to=max_val, by=increment)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*intxn_2
  
  # Compute variances
  var_1 = covMat[rhs1,rhs1] + (intxn_2^2)*covMat[paste0(rhs1, ':', intxn), paste0(rhs1, ':', intxn)] + 2*intxn_2*covMat[rhs1,  paste0(rhs1, ':', intxn)]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  d <- as.data.table(cbind(intxn_2, delta_1, upper_bound, lower_bound))
  int <- data.frame(inter = mod_frame[[intxn]])
  
  #Plot
  ggplot(environment = environment()) +
    ggtitle(title) +
    geom_line(data=d, aes(y=delta_1, x=intxn_2), lwd=1.5, color='black') +
    geom_line(data=d, aes(y=upper_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_line(data=d, aes(y=lower_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_hline(yintercept=0, lty=3) +
    #geom_rug(data=int, mapping = aes(x=inter), alpha=.05, color=secondcolor) +
    ylab(ylabel) +
    xlab(xlabel) +
    coord_cartesian(ylim = c(miny,maxy), xlim=c(min_val,max_val)) +
    theme(title = element_text(size=16, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.5),
          axis.title.x = element_text(size=18, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=18,angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=16,angle=0, vjust=.5, face="bold") ,
          axis.text.x = element_text(size=16, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(),
          legend.title=element_blank()
    )
  
}

drought <- continuousinteraction(felmmodel=drought.somefelm_fallow,
                                 lmmodel=drought.lm_fallow,
                                 rhs1="prcp.dev",
                                 intxn="temp.dev",
                                 minimum="min",
                                 maximum="max",
                                 miny=-0.1,
                                 maxy=0.1,
                                 incr="default",
                                 num_points = 1000,
                                 conf=.95,
                                 title="Drought (2006-2010)",
                                 xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                 ylabel="Effect of \u0394 Precipitation on Fallow Agriculture",
                                 maincolor='darkorange',
                                 secondcolor='darkgray')

predrought <- continuousinteraction(felmmodel=predrought.somefelm_fallow,
                                    lmmodel=predrought.lm_fallow,
                                    rhs1="prcp.dev",
                                    intxn="temp.dev",
                                    minimum="min",
                                    maximum="max",
                                    miny=-0.1,
                                    maxy=0.1,
                                    incr="default",
                                    num_points = 1000,
                                    conf=.95,
                                    title="Pre-Drought (2000-2005)",
                                    xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                    ylabel="Effect of \u0394 Precipitation on Fallow Agriculture",
                                    maincolor='darkorange',
                                    secondcolor='darkgray')

all <- continuousinteraction(felmmodel=all.somefelm_fallow,
                             lmmodel=all.lm_fallow,
                             rhs1="prcp.dev",
                             intxn="temp.dev",
                             minimum="min",
                             maximum="max",
                             miny=-0.1,
                             maxy=0.1,
                             incr="default",
                             num_points = 1000,
                             conf=.95,
                             title="Full Sample (2000-2010)",
                             xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                             ylabel="Effect of \u0394 Precipitation on Fallow Agriculture",
                             maincolor='darkorange',
                             secondcolor='darkgray')
## plot figure A1
cairo_pdf(file = 'climate_fallowness_intxn.pdf', width = 13, height = 6.5)
grid.arrange(all, predrought, drought, ncol=3)
dev.off()



###Appendix 2.1. Table A7: Testing Association between Drought error term Fallow on Migration Models###

all.fallow.resid <- model.frame(all.fallow.felm)
all.fallow.resid$all_fallow_resid <- all.fallow.felm$residuals

names(all.fallow.resid)[names(all.fallow.resid) == "name_en"] <- "name"
names(all.fallow.resid)[names(all.fallow.resid) == "as.factor(year)"] <- "year"
names(all.fallow.resid)[names(all.fallow.resid) == "as.factor(admin1)"] <- "admin1"
all.fallow.resid = subset(all.fallow.resid, select = -admin1 )
all.fallow.resid = subset(all.fallow.resid, select = -admin2 )
all.fallow.resid = subset(all.fallow.resid, select = -lights )
all.fallow.resid = subset(all.fallow.resid, select = -fallow )

all.fallow.resid$year = as.numeric(as.character(all.fallow.resid$year))
all.fallow.resid$all_fallow_resid = as.numeric(all.fallow.resid$all_fallow_resid)

master <- master %>%
  left_join(all.fallow.resid)

drought.fallow.resid <- model.frame(drought.fallow.felm)
drought.fallow.resid$drought_fallow_resid <- drought.fallow.felm$residuals

names(drought.fallow.resid)[names(drought.fallow.resid) == "name_en"] <- "name"
names(drought.fallow.resid)[names(drought.fallow.resid) == "as.factor(year)"] <- "year"
names(drought.fallow.resid)[names(drought.fallow.resid) == "as.factor(admin1)"] <- "admin1"
drought.fallow.resid = subset(drought.fallow.resid, select = -admin1 )
drought.fallow.resid = subset(drought.fallow.resid, select = -admin2 )
drought.fallow.resid = subset(drought.fallow.resid, select = -lights )
drought.fallow.resid = subset(drought.fallow.resid, select = -fallow )

drought.fallow.resid$year = as.numeric(as.character(drought.fallow.resid$year))
drought.fallow.resid$drought_fallow_resid = as.numeric(drought.fallow.resid$drought_fallow_resid)

master <- master %>%
  left_join(drought.fallow.resid)

predrought.fallow.resid <- model.frame(predrought.fallow.felm)
predrought.fallow.resid$predrought_fallow_resid <- predrought.fallow.felm$residuals

names(predrought.fallow.resid)[names(predrought.fallow.resid) == "name_en"] <- "name"
names(predrought.fallow.resid)[names(predrought.fallow.resid) == "as.factor(year)"] <- "year"
names(predrought.fallow.resid)[names(predrought.fallow.resid) == "as.factor(admin1)"] <- "admin1"
predrought.fallow.resid = subset(predrought.fallow.resid, select = -admin1 )
predrought.fallow.resid = subset(predrought.fallow.resid, select = -admin2 )
predrought.fallow.resid = subset(predrought.fallow.resid, select = -lights )
predrought.fallow.resid = subset(predrought.fallow.resid, select = -fallow )

predrought.fallow.resid$year = as.numeric(as.character(predrought.fallow.resid$year))
predrought.fallow.resid$predrought_fallow_resid = as.numeric(predrought.fallow.resid$predrought_fallow_resid)

master <- master %>%
  left_join(predrought.fallow.resid)

all.lm_error <- felm(all_fallow_resid ~ prcp.dev + temp.dev + prcp.dev:temp.dev| 0 | 0 | admin2, data=master)
drought.lm_error <- felm(drought_fallow_resid ~ prcp.dev + temp.dev + prcp.dev:temp.dev | 0 | 0 | admin2, data=master[master$year >= 2006 & master$year <= 2010,])
predrought.lm_error <- felm(predrought_fallow_resid ~ prcp.dev + temp.dev + prcp.dev:temp.dev | 0 | 0 | admin2, data=master[master$year < 2008,])


summary(all.lm_error)
summary(drought.lm_error)
summary(predrought.lm_error)

stargazer(all.lm_error, predrought.lm_error, drought.lm_error,
          type='latex',
          title = 'Climatic Stress and Residuals from Naive Fallowness Models',
          model.numbers = TRUE,
          dep.var.labels.include = FALSE,
          column.labels = c('All Years', 'Pre-Drought', 'Drought'),
          dep.var.caption = 'DV: Residuals from Time-by-Location FE Model of Fallowness on Nighttime Lights',
          covariate.labels = c('$\\Delta$ Precip.', '$\\Delta$ Temp.', '$\\Delta$ Precip. $*$ $\\Delta$ Temp.'),
          add.lines = list(c("Third-Level Admin FE", "No", "No", "No"),
                           c("Year FE", "No", "No", "No"),
                           c("First-Level Admin:Year FE", "No", "No", "No")),
          notes = c('Only the 185 nawahi Eklund et. al. 2022', 'record as having agricultural land',
                    'are included in the model.'),
          df=FALSE,
          header=FALSE,
          digits=3,
          no.space=T,
          font.size=NULL,
          column.sep.width="40pt"
)

######Figure A2: Visualizing Interaction between temperature and precipitation on residual errors from naive model #######

conf=.95
title="OLS of Climatic Stress on Residual Errors from Naive Model (2008-10)"
xlabel=expression(Delta ~ "Temperature (t-1)")
ylabel=expression("Effect of" ~Delta~" Precipitation (t-1)")

# Get coefficients of variables
beta_1 = -0.08600
beta_3 =  0.09679

# Set range of the moderator variable
# Minimum

min_val = 0
max_val = 2



# Create list of moderator values at which marginal effect is evaluated
x_2 <- seq(from=min_val, to=max_val, by=1)

# Compute marginal effects
delta_1 = beta_1 + beta_3*x_2

# Compute variances
var_1 =       0.001597363 + (x_2^2)*   0.001778621  + 2*x_2*   0.001778621    

# Standard errors
se_1 = sqrt(var_1)

# Upper and lower confidence bounds
z_score = qnorm(1 - ((1 - conf)/2))
upper_bound = delta_1 + z_score*se_1
lower_bound = delta_1 - z_score*se_1

# Determine the bounds of the graphing area
max_y = max(upper_bound)
min_y = min(lower_bound)

# Initialize plotting window

axis(side = 1, 
     at = 0:2,
     labels=c("0","1","2"))

plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title, xaxt='n')


# Plot estimated effects
lines(y=delta_1, x=x_2)
lines(y=upper_bound, x=x_2, lty=2)
lines(y=lower_bound, x=x_2, lty=2)

# Add a dashed horizontal line for zero
abline(h=0, lty=3)








###Figure A3####

continuousinteraction <- function(felmmodel, lmmodel, rhs1, intxn, minimum="min", maximum="max", miny=-8, maxy=1, incr="default", num_points = 1000, conf=.95, title="Interaction marginal effects plot", xlabel="Value of Intxn", ylabel="Estimated Marginal Effect", maincolor='black', secondcolor='black'){
  
  # Get coefficients
  beta_1 = felmmodel$coefficients[row.names(felmmodel$coefficients)==rhs1]
  
  # Get coefficient if interaction term==100
  beta_3 = felmmodel$coefficients[row.names(felmmodel$coefficients)==paste0(rhs1, ':', intxn)]
  
  # Extract Variance Covariance matrix
  covMat = felmmodel$clustervcv
  
  # Extract the data frame of the model
  mod_frame = model.frame(lmmodel)
  
  # Set range of the interaction variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[intxn]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[intxn]])
  }else{
    max_val = maximum
  }
  
  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum X value greater than maximum value.")
  }
  
  # Determine intervals between values of the x variable
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }
  
  # Create list of intxn values at which marginal effect is evaluated
  intxn_2 <- seq(from=min_val, to=max_val, by=increment)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*intxn_2
  
  # Compute variances
  var_1 = covMat[rhs1,rhs1] + (intxn_2^2)*covMat[paste0(rhs1, ':', intxn), paste0(rhs1, ':', intxn)] + 2*intxn_2*covMat[rhs1,  paste0(rhs1, ':', intxn)]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  d <- as.data.table(cbind(intxn_2, delta_1, upper_bound, lower_bound))
  int <- data.frame(inter = mod_frame[[intxn]])
  
  #Plot
  ggplot(environment = environment()) +
    ggtitle(title) +
    geom_line(data=d, aes(y=delta_1, x=intxn_2), lwd=1.5, color='black') +
    geom_line(data=d, aes(y=upper_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_line(data=d, aes(y=lower_bound, x=intxn_2), lty=2, lwd=1.5, color=maincolor) +
    geom_hline(yintercept=0, lty=3) +
    #geom_rug(data=int, mapping = aes(x=inter), alpha=.05, color=secondcolor) +
    ylab(ylabel) +
    xlab(xlabel) +
    coord_cartesian(ylim = c(miny,maxy), xlim=c(min_val,max_val)) +
    theme(title = element_text(size=16, angle=0, vjust=.5, face="bold"),
          plot.title = element_text(hjust=0.5),
          axis.title.x = element_text(size=18, angle=0, vjust=-0, face="bold"),
          axis.title.y = element_text(size=18,angle=90, vjust=1, face="bold"),
          axis.text.y = element_text(size=16,angle=0, vjust=.5, face="bold") ,
          axis.text.x = element_text(size=16, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line = element_line(),
          legend.title=element_blank()
    )
  
}

drought <- continuousinteraction(felmmodel=drought.felm_full,
                                 lmmodel=drought.lm_full,
                                 rhs1="prcp.dev",
                                 intxn="temp.dev",
                                 minimum="min",
                                 maximum="max",
                                 miny=-0.2,
                                 maxy=0.4,
                                 incr="default",
                                 num_points = 1000,
                                 conf=.95,
                                 title="Drought (2006-2010)",
                                 xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                 ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                                 maincolor='darkorange',
                                 secondcolor='darkgray')

predrought <- continuousinteraction(felmmodel=predrought.felm_full,
                                    lmmodel=predrought.lm_full,
                                    rhs1="prcp.dev",
                                    intxn="temp.dev",
                                    minimum="min",
                                    maximum="max",
                                    miny=-0.2,
                                    maxy=0.4,
                                    incr="default",
                                    num_points = 1000,
                                    conf=.95,
                                    title="Pre-Drought (2000-2005)",
                                    xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                                    ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                                    maincolor='darkorange',
                                    secondcolor='darkgray')

all <- continuousinteraction(felmmodel=all.felm_full,
                             lmmodel=all.lm_full,
                             rhs1="prcp.dev",
                             intxn="temp.dev",
                             minimum="min",
                             maximum="max",
                             miny=-0.2,
                             maxy=0.4,
                             incr="default",
                             num_points = 1000,
                             conf=.95,
                             title="Full Sample (2000-2010)",
                             xlabel="Mean Temperature \u0394 (in \u00B0 C)",
                             ylabel="Effect of \u0394 Precipitation on Nighttime Lights",
                             maincolor='darkorange',
                             secondcolor='darkgray')
## plot figure A3
cairo_pdf(file = 'climate_lights_intxn_fallow.pdf', width = 13, height = 6.5)
grid.arrange(all, predrought, drought, ncol=3)
dev.off()

